home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / Mathematics / Notebooks / SigProc2.0 / Packages / SignalProcessing / Digital / DSPAnalyze.m < prev    next >
Encoding:
Text File  |  1992-08-18  |  11.8 KB  |  393 lines
  1. (*  :Title:    Digital Signal Processing Analysis  *)
  2.  
  3. (*  :Authors:    Brian Evans, James McClellan  *)
  4.  
  5. (*
  6.     :Summary:    Perform basic analysis on a one-dimensional discrete-time
  7.         function, including plotting f[n] and graphing the poles
  8.         and zeroes of F[z].
  9.  *)
  10.  
  11. (*  :Context:    SignalProcessing`Digital`DSPAnalyze`  *)
  12.  
  13. (*  :PackageVersion:  2.7    *)
  14.  
  15. (*
  16.     :Copyright:    Copyright 1989-1991 by Brian L. Evans
  17.         Georgia Tech Research Corporation
  18.  
  19.     Permission to use, copy, modify, and distribute this software
  20.     and its documentation for any purpose and without fee is
  21.     hereby granted, provided that the above copyright notice
  22.     appear in all copies and that both that copyright notice and
  23.     this permission notice appear in supporting documentation,
  24.     and that the name of the Georgia Tech Research Corporation,
  25.     Georgia Tech, or Georgia Institute of Technology not be used
  26.     in advertising or publicity pertaining to distribution of the
  27.     software without specific, written prior permission.  Georgia
  28.     Tech makes no representations about the suitability of this
  29.     software for any purpose.  It is provided "as is" without
  30.     express or implied warranty.
  31.  *)
  32.  
  33. (*  :History:    *)
  34.  
  35. (*  :Keywords:    analysis, transforms, signals and systems    *)
  36.  
  37. (*  :Source:    *)
  38.  
  39. (*  :Warning:    *)
  40.  
  41. (*  :Mathematica Version:  1.2 or 2.0  *)
  42.  
  43. (*  :Limitation:  *)
  44.  
  45. (*  :Discussion:  *)
  46.  
  47. (*  :Functions:    DSPAnalyze  *)
  48.  
  49.  
  50.  
  51. (*  B E G I N     P A C K A G E  *)
  52.  
  53. BeginPackage[ "SignalProcessing`Digital`DSPAnalyze`",
  54.           "SignalProcessing`Digital`DTFT`",
  55.           "SignalProcessing`Digital`DSupport`",
  56.           "SignalProcessing`Digital`ZTransform`",
  57.           "SignalProcessing`Digital`ZSupport`",
  58.           "SignalProcessing`Support`TransSupport`",
  59.           "SignalProcessing`Support`ROC`",
  60.           "SignalProcessing`Support`SigProc`",
  61.           "SignalProcessing`Support`SupCode`" ]
  62.  
  63.  
  64. If [ TrueQ[ $VersionNumber >= 2.0 ],
  65.      Off[ General::spell ];
  66.      Off[ General::spell1 ] ];
  67.  
  68.  
  69. (*  U S A G E     I N F O R M A T I O N  *)
  70.  
  71. DSPAnalyze::usage =
  72.     "DSPAnalyze[f, n] and DSPAnalyze[f, n, start, end] \
  73.     will plot the real values of f[n] from n = Floor[start] to \
  74.     n = Ceiling[end], print its z-transform, indicate conditions \
  75.     for stability, display the pole-zero diagram, and plot the magnitude \
  76.     and phase responses. \
  77.     DSPAnalyze[f, {n1,n2}] and DSPAnalyze[f, {n1,n2}, {s1,s2}, {e1,e2}] \
  78.     treat f as a function of two-variables n1 and n2. \
  79.     The domain of the time-domain graph is n1 from s1 to e1 and \
  80.     n2 from s2 to e2. \
  81.     All non-variable symbols will be assigned a value of 1 when the \
  82.     analyzer generates graphics (override this with options like a -> 2)."
  83.  
  84. (*  E N D     U S A G E     I N F O R M A T I O N  *)
  85.  
  86.  
  87. Begin["`Private`"]
  88.  
  89.  
  90. (*  M E S S A G E S  *)
  91.  
  92. ZTransform::notvalid = "The forward z-transform could not be found."
  93. DSPAnalyze::hasval =
  94.     "The variable '' has value so it could not be used as a \
  95.     transform variable; '' was used instead."
  96. DTFTransform::notvalid = "The forward DTFT could not be found."
  97.  
  98.  
  99. (*  S U P P O R T I N G     F U N C T I O N S  *)
  100.  
  101. (*  dspAnalyzeQ  *)
  102. dspAnalyzeQ[start_, end_] :=
  103.     NumberQ[N[start]] && NumberQ[N[end]] && N[end >= start]
  104.  
  105. (*  checkvar  *)
  106. SetAttributes[checkvar, {HoldFirst}]
  107. checkvar[ xsym_, xval_, str_ ] := xsym /; SameQ[xsym, xval]
  108. checkvar[ xsym_, xval_, str_ ] :=
  109.     Block [ {sym},
  110.         sym = Unique["z"];
  111.         Message[ DSPAnalyze::hasval, str, sym ];
  112.         sym ]
  113.  
  114. (*  formatting-- this works if and only if the variables do not appear  *)
  115. (*      in a Block statement; i.e., they must be global to the package  *)
  116. w/:  Format[ w ] := "w"
  117. w1/: Format[ w1 ] := "w1"
  118. w2/: Format[ w2 ] := "w2"
  119.  
  120. z/:  Format[ z ] := "z"
  121. z1/: Format[ z1 ] := "z1"
  122. z2/: Format[ z2 ] := "z2"
  123.  
  124. (*  printDefaults  *)
  125. printDefaults[oplist_, varlist_, defaultvalues_] :=
  126.     Block [    {},
  127.         If [ EmptyQ[oplist],
  128.              Print[ "For plotting only, these symbols will" ];
  129.              Print[ "be assigned a value of 1:  ", varlist ],
  130.              Print[ "For plotting only, these symbols will" ];
  131.              Print[ "will be assigned default values:" ];
  132.              Print[ defaultvalues ] ] ]
  133.  
  134. (*  printStability  *)
  135. printStability[ztrans_] :=
  136.     Block [    {stable},
  137.         stable = Stable[ztrans];
  138.         Which [ SameQ[stable, True],
  139.               Print[ "The system is stable." ],
  140.             SameQ[stable, False],
  141.               Print[ "The system is not stable." ],
  142.             True,
  143.                Print[ "The system is stable if ", stable ] ];
  144.         stable ]
  145.  
  146. (*  zToFreqResp  *)
  147. zToFreqResp[ funct_, ztrans_, defaultvalues_, onedflag_, zrules_ ] :=
  148.     Block [ {mag, plotrange, rm, rp, stable, zfreqresp},
  149.  
  150.         (*  Print the sequence and its z-transform  *)
  151.         Print[ funct ];
  152.         Print[ "has the following z-transform:" ];
  153.         Print[ TheFunction[ztrans] ];
  154.  
  155.         (*  Print the Region of Convergence  *)
  156.         Print[ "The region of convergence is:" ];
  157.         rm = GetRMinus[ztrans];
  158.         rp = GetRPlus[ztrans];
  159.         If [ onedflag,
  160.              Print[ rm, " < |z| < ", rp ],
  161.              Print[ rm[[1]], " < |z1| < ", rp[[1]] ];
  162.              Print[ rm[[2]], " < |z2| < ", rp[[2]] ] ];
  163.  
  164.         (*  Print stability information  *)
  165.         stable = printStability[ztrans];
  166.  
  167.         (*  From now on, the default values are used  *)
  168.         If [ ! EmptyQ[defaultvalues],
  169.              Print[ "The default values are now being considered." ] ];
  170.  
  171.         (*  Show the Pole-Zero Plot  *)
  172.         PoleZeroPlot[ztrans /. defaultvalues];
  173.  
  174.         If [ stable /. defaultvalues,
  175.              Print[ "Since the sequence is stable, the" ];
  176.              Print[ "frequency response can be computed" ];
  177.              Print[ "directly from the z-transform:" ];
  178.              zfreqresp = TheFunction[ztrans] /. zrules;
  179.              Print[ zfreqresp ];
  180.              plotrange = All;
  181.              mag = Log,
  182.  
  183.              Print[ "Since the sequence is not stable, the" ];
  184.              Print[ "frequency response has no meaning, but" ];
  185.              Print[ "it will be plotted anyway.  However, the DC" ];
  186.              Print[ "(zero frequency) term will be Infinity which" ];
  187.              Print[ "will generate a warning from Mathematica." ];
  188.              plotrange = Automatic;
  189.              mag = Linear,
  190.               
  191.              Print[ "Even though the stability of the sequence" ];
  192.              Print[ "is not known, the frequency response will" ];
  193.              Print[ "still be plotted." ];
  194.              plotrange = Automatic;
  195.              mag = Linear ];
  196.  
  197.         { mag, plotrange, zfreqresp } ]
  198.  
  199.  
  200. (*  D S P A n a l y z e  *)
  201.  
  202. DSPAnalyze[f_, n_Symbol] := DSPAnalyze[f, n, 0]
  203. DSPAnalyze[f_, n_Symbol, start_] := DSPAnalyze[f, n, start, start + 40]
  204.  
  205. DSPAnalyze[f_, n_Symbol, start_, end_, op___] :=
  206.     Block [    {bogus, defaultvalues = {}, dtft, freqresp, fnumeric,
  207.          ftemp, funct, mag, oplist, plotrange = All, varlist,
  208.          zfreqresp, ztrans},
  209.  
  210.         oplist = ToList[op];
  211.         funct = ToDiscrete[f];
  212.  
  213.         (*  DISCRETE-TIME ANALYSIS *)
  214.  
  215.                 (*  Convert expression to plottable/functional form    *)
  216.         (*  and replace constants like Pi with numerical value *)
  217.  
  218.         fnumeric = N[ TheFunction[ funct ] ];
  219.  
  220.                 (*  Find all valueless symbols besides the time variable *)
  221.         (*  The first parameter in Summation is a local symbol   *)
  222.  
  223.         ftemp = fnumeric /.
  224.               ( Summation[m_, l_, u_, inc_][expr_] :>
  225.                 bogus[l, u, inc, GetVariables[expr /. m -> l]] );
  226.         varlist = Select[GetVariables[ftemp /. oplist],
  227.                  Function[var, ! SameQ[var, n]]];
  228.         If [ (! EmptyQ[varlist]) || (! EmptyQ[oplist]),
  229.              defaultvalues = oplist ~Join~ Map[(#1 -> 1)&, varlist];
  230.              printDefaults[oplist, varlist, defaultvalues];
  231.              fnumeric = fnumeric /. defaultvalues ];
  232.  
  233.         (*  Plot "time"-domain    *)
  234.  
  235.         SequencePlot[ fnumeric, {n, start, end},
  236.                   PlotLabel -> "Discrete-Time Domain Analysis" ];
  237.  
  238.         (*  Z-DOMAIN ANALYSIS *)
  239.  
  240.         (*    Find F[z]; if it exists, print it out, give stability  *)
  241.         (*  information, and display a pole-zero plot.             *)
  242.  
  243.         ztrans = ZTransform[funct, n, z];
  244.         If [ ! SameQ[Head[ztrans], ZTransData],
  245.              Message[ZTransform::notvalid];
  246.              mag = Linear,
  247.              { mag, plotrange, zfreqresp } =
  248.                zToFreqResp[ funct, ztrans, defaultvalues,
  249.                     True, (z -> Exp[I w]) ] ];
  250.  
  251.         (*  FREQUENCY-DOMAIN ANALYSIS  *)
  252.  
  253.         (*  Determine the frequency response  *)
  254.  
  255.         dtft = DTFTransform[funct, n, w, Simplify -> True];
  256.         If [ ! SameQ[Head[dtft], DTFTData],
  257.              Message[DTFTransform::notvalid]; Return[ztrans] ];
  258.         freqresp = TheFunction[dtft];
  259.  
  260.         (*  Plot the interesting section of the freq. response  *)
  261.  
  262.         If [ ! SameQ[freqresp, zfreqresp],
  263.              Print[ funct ];
  264.              Print[ "has the following frequency response:"];
  265.              Print[ freqresp ];
  266.              Print[ "where the frequency response is periodic" ];
  267.              Print[ "usually with period 2 Pi" ] ];
  268.  
  269.         freqresp = ToContinuous[freqresp];
  270.         MagPhasePlot[ freqresp /. defaultvalues,
  271.                   {w, -2 Pi, 2 Pi},
  272.                   Domain -> Continuous, DomainScale -> Linear,
  273.                   MagRangeScale -> mag, PhaseRangeScale -> Degree,
  274.                   PlotRange -> plotrange ];
  275.  
  276.         (*  Return the z-transform  *)
  277.  
  278.         ztrans /. ( z -> checkvar[Global`z, Global`z, "z"] ) ] /;
  279.     N [ dspAnalyzeQ[start, end] ]
  280.  
  281.  
  282. DSPAnalyze[f_, {n1_Symbol, n2_Symbol}] := DSPAnalyze[f, {n1, n2}, {0, 0}]
  283. DSPAnalyze[f_, {n1_Symbol, n2_Symbol}, start_] :=
  284.     DSPAnalyze[f, {n1, n2}, start, start + {9, 9}]
  285.  
  286. DSPAnalyze[f_, {n1_Symbol, n2_Symbol}, start_, end_, op___] :=
  287.     Block [    {coordlist, defaultvalues = {}, dtft, freqresp, ftemp, funct,
  288.          mag, nlist, oplist, plotrange = All, stable, varlist,
  289.          wlist, zfreqresp, zlist, ztrans},
  290.  
  291.         oplist = ToList[op];
  292.  
  293.         nlist = { n1, n2 };
  294.         wlist = { w1, w2 };
  295.         zlist = { z1, z2 };
  296.         funct = ToDiscrete[f];
  297.  
  298.         (*  DISCRETE-TIME ANALYSIS *)
  299.  
  300.                 (*  Convert expression to plottable/functional form    *)
  301.         (*  and replace constants like Pi with numerical value *)
  302.  
  303.         fnumeric = N[ TheFunction[ funct ] ];
  304.  
  305.                 (*  Find all valueless symbols besides the time variable *)
  306.         (*  The first parameter in Summation is a local symbol   *)
  307.  
  308.         ftemp = fnumeric /.
  309.               ( Summation[m_, l_, u_, inc_][expr_] :>
  310.                 bogus[l, u, inc, GetVariables[expr /. m -> l]] );
  311.         varlist = Select[GetVariables[ftemp /. oplist],
  312.                   Function[var,
  313.                       ! SameQ[var,n1] && ! SameQ[var,n2]]];
  314.         If [ (! EmptyQ[varlist]) || (! EmptyQ[oplist]),
  315.              defaultvalues = oplist ~Join~ Map[(#1 -> 1)&, varlist];
  316.              printDefaults[oplist, varlist, defaultvalues];
  317.              fnumeric = fnumeric /. defaultvalues ];
  318.  
  319.         (*  Generate three-dimensional coordinates for plot  *)
  320.  
  321.         SequencePlot[ fnumeric,
  322.                   { n1, start[[1]], end[[1]] },
  323.                   { n2, start[[2]], end[[2]] },
  324.                   PlotLabel -> "Discrete-Time Domain Analysis" ];
  325.  
  326.         (*  Z-DOMAIN ANALYSIS  *)
  327.  
  328.         (*  Find F[z] and print it out.  *)
  329.  
  330.         ztrans = ZTransform[funct, nlist, zlist];
  331.         If [ ! SameQ[Head[ztrans], ZTransData],
  332.              Message[ZTransform::notvalid];
  333.              mag = Linear,
  334.  
  335.              { mag, plotrange, zfreqresp } =
  336.                zToFreqResp[ funct, ztrans, defaultvalues, False,
  337.                 ReplaceWith[zlist, {Exp[I w1], Exp[I w2]}] ] ];
  338.  
  339.         (*  FREQUENCY-DOMAIN ANALYSIS  *)
  340.  
  341.         (*  Determine the frequency response  *)
  342.  
  343.         dtft = DTFTransform[funct, nlist, wlist];
  344.         If [ ! SameQ[Head[dtft], DTFTData],
  345.              Message[DTFT::notvalid]; Return[ztrans] ];
  346.         freqresp = TheFunction[dtft];
  347.  
  348.         (*  Plot the interesting section of the freq. response  *)
  349.  
  350.         If [ ! SameQ[freqresp, zfreqresp],
  351.              Print[ funct ];
  352.              Print[ "has the following frequency response:" ];
  353.              Print[ freqresp ] ];
  354.  
  355.         (*  Making the freqresp periodic takes too long to plot: *)
  356.         (*    Periodic[{2 Pi, 2 Pi}, {w1, w2}][freqresp]     *)
  357.  
  358.         MagPhasePlot[ freqresp /. defaultvalues,
  359.                   {w1, 0, 2 Pi}, {w2, 0, 2 Pi},
  360.                   Domain -> Continuous, DomainScale -> Linear,
  361.                   MagRangeScale -> mag, PhaseRangeScale -> Degree,
  362.                   PlotRange -> plotrange ];
  363.  
  364.         (*  Return the z-transform  *)
  365.  
  366.         ztrans /. { z1 -> checkvar[Global`z1, Global`z1, "z1"],
  367.                 z2 -> checkvar[Global`z2, Global`z2, "z2"] } ] /;
  368.  
  369.     N [ dspAnalyzeQ[ start[[1]], end[[1]] ] &&
  370.         dspAnalyzeQ[ start[[2]], end[[2]] ] ]
  371.  
  372.  
  373. (*  E N D     P A C K A G E  *)
  374.  
  375. End[]
  376. EndPackage[]
  377.  
  378. If [ TrueQ[ $VersionNumber >= 2.0 ],
  379.      On[ General::spell ];
  380.      On[ General::spell1 ] ];
  381.  
  382.  
  383. (*  H E L P     I N F O R M A T I O N  *)
  384.  
  385. Combine[ SPfunctions, { DSPAnalyze } ]
  386. Protect[ DSPAnalyze ]
  387.  
  388.  
  389. (*  E N D I N G     M E S S A G E  *)
  390.  
  391. Print["The digital signal analyzer is now loaded."]
  392. Null
  393.